require(runjags)
require(coda)
require(stringi)

source("A10_functions.r")

data <- read.csv("C1_main_data.csv")
data$date <- as.Date(data$date)

#### estimation ####
# preprocessing
n <- nrow(data) - 1
m <- 4
y <- cbind(data$K_L_3s, data$K_D_3s, data$K_C_3s, data$K_K_3s)[-1, ]
k <- rowSums(y)
y[k == 0, ] <- NA
k[k == 0] <- 1

x <- cbind(diff(data$cabinet / 100), 
           diff(data$LDP / 100), 
           data$DPJ.gov[-1], 
           diff(data$cabinet / 100) * data$DPJ.gov[-1], 
           diff(data$LDP / 100) * data$DPJ.gov[-1], 
           data$HOR[-1] / 47, 
           (data$HOR[-1] / 47) ^ 2, 
           data$HOC[-1] / 35, 
           (data$HOC[-1] / 35) ^ 2, 
           diff(data$pmid))
l <- ncol(x)

LDP.mean.minus.SD <- mean(x[, 2]) - sd(x[, 2])
LDP.mean.plus.SD <- mean(x[, 2]) + sd(x[, 2])

# initial values
initial.values <- list()
for (i in 1:3) {
  initial.values[[i]] <- inits.fun(10 * i, 100 * i, i)
}

# estimation
result <- run.jags(model = "B1_dynamic_multinomial_model.R", 
                   monitor = c("beta", "Omega", "theta", "theta.star"), 
                   data = list(y = y, k = k, x = x, n = n, m = m, l = l, 
                               mu0 = rep(0, m - 1), 
                               Omega.star0 = diag(m - 1) * 0.01), 
                   inits = initial.values, modules = "glm", 
                   n.chains = 3, burnin = 5000, sample = 5000, 
                   adapt = 20000, summarise = FALSE, 
                   thin = 50, method = "parallel")
# save(result, file = "D4_Komeito_interaction_result.Rdata")
# load("D4_Komeito_interaction_result.Rdata")

# convergence check
sum(gelman.diag(result$mcmc, multivariate = FALSE)$psrf[, 1] > 1.1)

# combine MCMC draws into a single matrix
sims.matrix <- rbind(result$mcmc[[1]], result$mcmc[[2]], result$mcmc[[3]])
colnames(sims.matrix) <- colnames(result$mcmc[[1]])

#### results ####
var.names <- c("Lagged cabinet approval", "Lagged LDP support", 
               "DPJ government", "Cabinet x DPJ gov.", "LDP x DPJ gov.", 
               "Time since HoR election", 
               "Time since HoR election squared", "Time since HoC election", 
               "Time since HoC election squared", "Inauguration of the PM")

# Table A.6
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], 
                        var.names), 3)
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], 
                        var.names), 3)
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], 
                        var.names), 3)

#### simulation for the LDP support rate ####
LDP.seq <- seq(round(min(x[, 2]), 3), round(max(x[, 2]), 3), 0.001)

# Figure A.3
sim.result.LDP.gov <- array(NA, c(length(LDP.seq), 4, 3))
for (i in 1:length(LDP.seq)) {
  sim.x <- x
  sim.x[, 2] <- LDP.seq[i]
  sim.x[, 3] <- sim.x[, 4] <- sim.x[, 5] <- 0
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1 <- exp.theta.star.1 / sum.exp.theta.star
  prob.2 <- exp.theta.star.2 / sum.exp.theta.star
  prob.3 <- exp.theta.star.3 / sum.exp.theta.star
  prob.4 <- 1 / sum.exp.theta.star
  sim.result.LDP.gov[i, 1, ] <- posterior.summary(rowMeans(prob.1))
  sim.result.LDP.gov[i, 2, ] <- posterior.summary(rowMeans(prob.2))
  sim.result.LDP.gov[i, 3, ] <- posterior.summary(rowMeans(prob.3))
  sim.result.LDP.gov[i, 4, ] <- posterior.summary(rowMeans(prob.4))
}

sim.result.LDP.opp <- array(NA, c(length(LDP.seq), 4, 3))
for (i in 1:length(LDP.seq)) {
  sim.x <- x
  sim.x[, 2] <- sim.x[, 5] <- LDP.seq[i]
  sim.x[, 3] <- 1
  sim.x[, 4] <- sim.x[, 1]
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1 <- exp.theta.star.1 / sum.exp.theta.star
  prob.2 <- exp.theta.star.2 / sum.exp.theta.star
  prob.3 <- exp.theta.star.3 / sum.exp.theta.star
  prob.4 <- 1 / sum.exp.theta.star
  sim.result.LDP.opp[i, 1, ] <- posterior.summary(rowMeans(prob.1))
  sim.result.LDP.opp[i, 2, ] <- posterior.summary(rowMeans(prob.2))
  sim.result.LDP.opp[i, 3, ] <- posterior.summary(rowMeans(prob.3))
  sim.result.LDP.opp[i, 4, ] <- posterior.summary(rowMeans(prob.4))
}

png("Figure_A3.png", 
    width = 6, height = 1.6, units = "in", pointsize = 7, res = 1500)
layout(matrix(1:4, 1, 4))
par(mar = c(4, 4, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0.05, 0.25), 
     main = "LDP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
line.fun.2(LDP.seq, sim.result.LDP.opp[, 1, ])
line.fun(LDP.seq, sim.result.LDP.gov[, 1, ])
rug(x[x[, 3] == 0, 2])
rug(x[x[, 3] == 1, 2], side = 3, col = "gray")
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0, 0.2), 
     main = "DPJ", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
line.fun.2(LDP.seq, sim.result.LDP.opp[, 2, ])
line.fun(LDP.seq, sim.result.LDP.gov[, 2, ])
rug(x[x[, 3] == 0, 2])
rug(x[x[, 3] == 1, 2], side = 3, col = "gray")
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0, 0.2), 
     main = "JCP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
line.fun.2(LDP.seq, sim.result.LDP.opp[, 3, ])
line.fun(LDP.seq, sim.result.LDP.gov[, 3, ])
rug(x[x[, 3] == 0, 2])
rug(x[x[, 3] == 1, 2], side = 3, col = "gray")
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0.4, 0.8), 
     main = "Komeito", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
line.fun.2(LDP.seq, sim.result.LDP.opp[, 4, ])
line.fun(LDP.seq, sim.result.LDP.gov[, 4, ])
rug(x[x[, 3] == 0, 2])
rug(x[x[, 3] == 1, 2], side = 3, col = "gray")
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
dev.off()